1 Colors plot

colors <- c("LMER" = "#DCA237", "Anova" = "#459B76")
lines <- c("Balanced" = "solid", "Unbalanced" = "dotted")

2 BALANCED DESIGN

3 Normal distribution

3.1 Simulation

# Power Analyses
nb_simul <- 10000

#Perform siluations
sim <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", 
              rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = 5, nfruit = 3, nhab = 3, nrep = 10, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim <- as.data.frame(t(sim))

#Add indic SA
sim <- add_indic_sign(sim)

#Add number of simul for each SA value to keep only 100 simuls per SA value
sim <- add_sim_number_SA(sim)


## Local adaptation (genetic)
tapply(sim$IndicGen[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], length)
## 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 1.3 1.4 1.5 
##  10 100 100 100 100 100 100 100 100 100 100 100  37  19   1
tapply(sim$IndicGen_aov[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], mean)
##       0.1       0.2       0.3       0.4       0.5       0.6       0.7       0.8 
## 0.6000000 0.5800000 0.6200000 0.5200000 0.6600000 0.6400000 0.6700000 0.7000000 
##       0.9         1       1.1       1.2       1.3       1.4       1.5 
## 0.6400000 0.7800000 0.7600000 0.8300000 0.8648649 0.8421053 1.0000000
## Local adaptation (non genetic)
tapply(sim$IndicNonGen[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##       0.1       0.2       0.3       0.4       0.5       0.6       0.7       0.8 
## 0.0000000 0.1100000 0.1800000 0.2000000 0.2300000 0.3100000 0.2800000 0.5400000 
##       0.9         1       1.1       1.2       1.3       1.4       1.5       1.7 
## 0.4900000 0.5400000 0.6600000 0.6938776 0.6346154 0.9230769 0.8000000 1.0000000
tapply(sim$IndicNonGen_aov[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##       0.1       0.2       0.3       0.4       0.5       0.6       0.7       0.8 
## 0.0000000 0.1400000 0.1800000 0.1900000 0.2300000 0.2900000 0.2700000 0.5400000 
##       0.9         1       1.1       1.2       1.3       1.4       1.5       1.7 
## 0.4900000 0.5300000 0.6600000 0.7040816 0.6346154 0.9230769 1.0000000 1.0000000
tapply(sim$IndicGen, sim$SA_Gen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5 
##   10  107  496 1127 1633 1919 1672 1307  868  464  234  106   37   19    1
tapply(sim$IndicNonGen, sim$SA_NonGen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.7 
##    4  111  493 1054 1759 1912 1695 1232  853  451  267   98   52   13    5    1

3.2 Plot

####################################
###### Plot genetic
####################################


#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5 
##   10  107  496 1127 1633 1919 1672 1307  868  464  234  106   37   19    1
#Save the same number of simuls for each SA value
nb_per_SA=100
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")
#To check 
# dim(sim)
# dim(sim_select_genetic)
# tapply(sim_select_genetic$index_SA_Gen, sim_select_genetic$SA_Gen, length)

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))



#Plot
plot_genetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Normal data") +
   theme_LO_sober  + 
   ggtitle("Normal data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_normal

# #Save plot
# name_plot<-paste0("Simul_powertest_genetic_normal_",
#       nb_simul,"nbsimul_",nb_per_SA,
#       "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   plot_genetic_normal, base_height = 10/cm(1), base_width = 14/cm(1), dpi = 1200)


################################# 
### Plot non-genetic
#################################
#Number of simuls for each SA non genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.7 
##    4  111  493 1054 1759 1912 1695 1232  853  451  267   98   52   13    5    1
#Save the same number of simuls for each SA value
nb_per_SA 
## [1] 100
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")
#To check 
# dim(sim)
# dim(sim_select_nongenetic)
# tapply(sim_select_nongenetic$index_SA_Gen, sim_select_nongenetic$SA_Gen, length)
# tapply(sim_select_nongenetic$index_SA_NonGen, sim_select_nongenetic$SA_NonGen, length)

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))


colors<-c("LMER" = "#DCA237", "Anova" = "#459B76")
#Plot
plot_nongenetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Normal data") +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_normal

#Save plot
# name_plot<-paste0("Simul_powertest_nongenetic_normal_",
#       nb_simul,"nbsimul_",nb_per_SA,
#       "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   plot_nongenetic_normal, base_height = 10/cm(1), base_width = 14/cm(1), dpi = 1200)


#################################
### Common plot
#################################

power_plot_normal<-cowplot::plot_grid(plot_nongenetic_normal,plot_genetic_normal,
                          labels=c("A", "B"), 
                          label_size = 15,
                          ncol =1, nrow = 2,
                     #    hjust = 0, vjust = 1,
                          scale = c(1,  1))
power_plot_normal

#Save plot
name_plot<-paste0("Simul_powertest_normal_",
      nb_simul,"nbsimul_",nb_per_SA,
      "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   power_plot_normal, base_height = 20/cm(1), base_width = 14/cm(1), dpi = 1200)

3.3 False positive rate analyses

### No genetic and no non genetic effects
nb_simul = 200
sim <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", rho=0, rho_ng=0, 
              npop_per_fruit = 5, nfruit = 3, nhab = 3, nrep = 10, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
  
## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.01
mean(sim$IndicGen_aov)
## [1] 0.01
## Local adaptation (non genetic)
mean(sim$IndicNonGen)
## [1] 0.02
mean(sim$IndicNonGen_aov)
## [1] 0.025
### Genetic effects and no non-genetic effects
sim <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", rho=-0.1, rho_ng=0,  
              npop_per_fruit = 5, nfruit = 3, nhab = 3, nrep = 10, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
tapply(sim$IndicGen, sim$SA_Gen, mean)
##        0.2        0.3        0.4        0.5        0.6        0.7        0.8 
## 0.00000000 0.16666667 0.06666667 0.35000000 0.15384615 0.13333333 0.30769231 
##        0.9          1        1.1        1.2        1.3 
## 0.26086957 0.54545455 0.66666667 1.00000000 1.00000000
tapply(sim$IndicGen_aov, sim$SA_Gen, mean)
##        0.2        0.3        0.4        0.5        0.6        0.7        0.8 
## 0.00000000 0.16666667 0.06666667 0.35000000 0.15384615 0.13333333 0.30769231 
##        0.9          1        1.1        1.2        1.3 
## 0.26086957 0.54545455 0.66666667 1.00000000 1.00000000
## Local adaptation (non genetic)
mean(sim$IndicNonGen)
## [1] 0.02
mean(sim$IndicNonGen_aov)
## [1] 0.025
### Non genetic effects and no genetic effects
sim <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", rho=0, rho_ng=-0.1,  
              npop_per_fruit = 5, nfruit = 3, nhab = 3, nrep = 10, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.08
mean(sim$IndicGen_aov)
## [1] 0.08
## Local adaptation (non genetic)
tapply(sim$IndicNonGen, sim$SA_NonGen, mean)
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.0000000 0.1250000 0.2000000 0.2187500 0.4285714 0.2702703 0.5000000 0.5333333 
##         1       1.1       1.2 
## 0.2500000 0.5000000 1.0000000
tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.0000000 0.2500000 0.1500000 0.2187500 0.4285714 0.2702703 0.5000000 0.5333333 
##         1       1.1       1.2 
## 0.2500000 0.5000000 1.0000000

4 Poisson distribution

4.1 Simulation

nb_simul <- 10000

#Perform siluations
sim <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",  rho = -0.05, rho_ng = -0.05, 
              npop_per_fruit = 5, nfruit = 3, nhab = 3, nrep = 10, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim <- as.data.frame(t(sim))

#Add indic SA
sim <- add_indic_sign(sim)

#Add number of simul for each SA value to keep only 100 simuls per SA value
sim <- add_sim_number_SA(sim)


## Local adaptation (genetic)
tapply(sim$IndicGen[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], length)
## 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 1.3 1.4 1.5 
##  10 100 100 100 100 100 100 100 100 100 100 100  37  19   1
tapply(sim$IndicGen_aov[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], mean)
##       0.1       0.2       0.3       0.4       0.5       0.6       0.7       0.8 
## 0.6000000 0.6400000 0.6600000 0.6100000 0.7200000 0.7000000 0.8300000 0.8400000 
##       0.9         1       1.1       1.2       1.3       1.4       1.5 
## 0.8300000 0.8900000 0.8900000 0.9400000 0.9189189 0.9473684 1.0000000
## Local adaptation (non genetic)
tapply(sim$IndicNonGen[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##       0.1       0.2       0.3       0.4       0.5       0.6       0.7       0.8 
## 0.0000000 0.0200000 0.1100000 0.1700000 0.3000000 0.4300000 0.4300000 0.6000000 
##       0.9         1       1.1       1.2       1.3       1.4       1.5       1.7 
## 0.6100000 0.6800000 0.8100000 0.8469388 0.7884615 0.8461538 0.8000000 1.0000000
tapply(sim$IndicNonGen_aov[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##       0.1       0.2       0.3       0.4       0.5       0.6       0.7       0.8 
## 0.2500000 0.1300000 0.2000000 0.2700000 0.3600000 0.4200000 0.4300000 0.5700000 
##       0.9         1       1.1       1.2       1.3       1.4       1.5       1.7 
## 0.6400000 0.6800000 0.8200000 0.8673469 0.7884615 0.8461538 0.8000000 1.0000000
tapply(sim$IndicGen, sim$SA_Gen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5 
##   10  107  496 1127 1633 1919 1672 1307  868  464  234  106   37   19    1
tapply(sim$IndicNonGen, sim$SA_NonGen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.7 
##    4  111  493 1054 1759 1912 1695 1232  853  451  267   98   52   13    5    1

4.2 Plot

####################################
###### Plot genetic
####################################


#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5 
##   10  107  496 1127 1633 1919 1672 1307  868  464  234  106   37   19    1
#Save the same number of simuls for each SA value
nb_per_SA=100
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")
#To check 
# dim(sim)
# dim(sim_select_genetic)
# tapply(sim_select_genetic$index_SA_Gen, sim_select_genetic$SA_Gen, length)

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))


#Plot
plot_genetic_poisson<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober  + 
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_poisson

# #Save plot
# name_plot<-paste0("Simul_powertest_genetic_poisson_",
#       nb_simul,"nbsimul_",nb_per_SA,
#       "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   plot_genetic_poisson, base_height = 10/cm(1), base_width = 14/cm(1), dpi = 1200)


################################# 
### Plot non-genetic
#################################
#Number of simuls for each SA non genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.7 
##    4  111  493 1054 1759 1912 1695 1232  853  451  267   98   52   13    5    1
#Save the same number of simuls for each SA value
nb_per_SA 
## [1] 100
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")
#To check 
# dim(sim)
# dim(sim_select_nongenetic)
# tapply(sim_select_nongenetic$index_SA_Gen, sim_select_nongenetic$SA_Gen, length)
# tapply(sim_select_nongenetic$index_SA_NonGen, sim_select_nongenetic$SA_NonGen, length)

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))


colors<-c("LMER" = "#DCA237", "Anova" = "#459B76")
#Plot
plot_nongenetic_poisson<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   ggtitle("Poisson data") +
   theme_LO_sober +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_poisson

#Save plot
# name_plot<-paste0("Simul_powertest_nongenetic_poisson_",
#       nb_simul,"nbsimul_",nb_per_SA,
#       "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   plot_nongenetic_poisson, base_height = 10/cm(1), base_width = 14/cm(1), dpi = 1200)


#################################
### Common plot
#################################

power_plot_poisson<-cowplot::plot_grid(plot_nongenetic_poisson,plot_genetic_poisson,
                          labels=c("A", "B"), 
                          label_size = 15,
                          ncol =1, nrow = 2,
                     #    hjust = 0, vjust = 1,
                          scale = c(1,  1))
power_plot_poisson

#Save plot
# name_plot<-paste0("Simul_powertest_poissondata_",
#       nb_simul,"nbsimul_",nb_per_SA,
#       "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   power_plot_poisson, base_height = 20/cm(1), base_width = 14/cm(1), dpi = 1200)
# 
# 

4.3 False positive rate analyses

### No genetic and no non genetic effects
nb_simul = 200
sim <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", rho=0, rho_ng=0,  
              npop_per_fruit = 5, nfruit = 3, nhab = 3, nrep = 10, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
  
## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.025
mean(sim$IndicGen_aov)
## [1] 0.025
## Local adaptation (non genetic)
mean(sim$IndicNonGen)
## [1] 0.01
mean(sim$IndicNonGen_aov)
## [1] 0.03
### Genetic effects and no non-genetic effects
sim <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", rho=-0.1, rho_ng=0,  
              npop_per_fruit = 5, nfruit = 3, nhab = 3, nrep = 10, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
tapply(sim$IndicGen, sim$SA_Gen, mean)
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.2500000 0.0000000 0.1333333 0.4500000 0.2564103 0.2333333 0.3846154 0.3478261 
##         1       1.1       1.2       1.3 
## 0.5454545 0.8333333 1.0000000 1.0000000
tapply(sim$IndicGen_aov, sim$SA_Gen, mean)
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.2500000 0.0000000 0.1333333 0.4500000 0.2564103 0.2333333 0.3846154 0.3478261 
##         1       1.1       1.2       1.3 
## 0.5454545 0.8333333 1.0000000 1.0000000
## Local adaptation (non genetic)
mean(sim$IndicNonGen)
## [1] 0.045
mean(sim$IndicNonGen_aov)
## [1] 0.05
### Non genetic effects and no genetic effects
sim <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", rho=0, rho_ng=-0.1, 
              npop_per_fruit = 5, nfruit = 3, nhab = 3, nrep = 10, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)

sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.095
mean(sim$IndicGen_aov)
## [1] 0.095
## Local adaptation (non genetic)
tapply(sim$IndicNonGen, sim$SA_NonGen, mean)
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.0000000 0.0000000 0.0500000 0.1250000 0.3673469 0.2432432 0.2916667 0.5333333 
##         1       1.1       1.2 
## 0.2500000 0.5000000 1.0000000
tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.0000000 0.0000000 0.1500000 0.1875000 0.2857143 0.2702703 0.3750000 0.4000000 
##         1       1.1       1.2 
## 0.2500000 0.5000000 1.0000000

5 Binomial distribution

5.1 Simulation

nb_simul <- 10000

#Perform siluations
sim <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",  rho = -0.05, 
              rho_ng = -0.05, 
              npop_per_fruit = 5, nfruit = 3, nhab = 3, nrep = 10, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5, ntrial = 60)
sim <- as.data.frame(t(sim))


#Add indic SA
sim <- add_indic_sign(sim)

#Add number of simul for each SA value to keep only 100 simuls per SA value
sim <- add_sim_number_SA(sim)


## Local adaptation (genetic)
tapply(sim$IndicGen[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], length)
## 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 1.3 1.4 1.5 
##  10 100 100 100 100 100 100 100 100 100 100 100  37  19   1
tapply(sim$IndicGen_aov[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], mean)
##       0.1       0.2       0.3       0.4       0.5       0.6       0.7       0.8 
## 0.6000000 0.6300000 0.6000000 0.5400000 0.6500000 0.6500000 0.7200000 0.7000000 
##       0.9         1       1.1       1.2       1.3       1.4       1.5 
## 0.6700000 0.8100000 0.8100000 0.8400000 0.8918919 0.8421053 1.0000000
## Local adaptation (non genetic)
tapply(sim$IndicNonGen[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##       0.1       0.2       0.3       0.4       0.5       0.6       0.7       0.8 
## 0.0000000 0.2000000 0.1600000 0.2200000 0.2700000 0.3500000 0.3700000 0.5300000 
##       0.9         1       1.1       1.2       1.3       1.4       1.5       1.7 
## 0.5500000 0.5400000 0.6500000 0.7040816 0.6538462 0.9230769 0.8000000 1.0000000
tapply(sim$IndicNonGen_aov[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)
##       0.1       0.2       0.3       0.4       0.5       0.6       0.7       0.8 
## 0.0000000 0.2000000 0.1700000 0.2100000 0.2700000 0.3500000 0.3700000 0.5400000 
##       0.9         1       1.1       1.2       1.3       1.4       1.5       1.7 
## 0.5500000 0.5400000 0.6500000 0.7040816 0.6538462 0.9230769 0.8000000 1.0000000
tapply(sim$IndicGen, sim$SA_Gen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5 
##   10  107  496 1127 1633 1919 1672 1307  868  464  234  106   37   19    1
tapply(sim$IndicNonGen, sim$SA_NonGen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.7 
##    4  111  493 1054 1759 1912 1695 1232  853  451  267   98   52   13    5    1

5.2 Plot

####################################
###### Plot genetic
####################################


#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5 
##   10  107  496 1127 1633 1919 1672 1307  868  464  234  106   37   19    1
#Save the same number of simuls for each SA value
nb_per_SA=100
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")
#To check 
# dim(sim)
# dim(sim_select_genetic)
# tapply(sim_select_genetic$index_SA_Gen, sim_select_genetic$SA_Gen, length)

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))


#Plot
plot_genetic_binomial<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Binomial data") +
   theme_LO_sober +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_binomial

# #Save plot
# name_plot<-paste0("Simul_powertest_genetic_binomial_",
#       nb_simul,"nbsimul_",nb_per_SA,
#       "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   plot_genetic_binomial, base_height = 10/cm(1), base_width = 14/cm(1), dpi = 1200)


################################# 
### Plot non-genetic
#################################
#Number of simuls for each SA non genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.7 
##    4  111  493 1054 1759 1912 1695 1232  853  451  267   98   52   13    5    1
#Save the same number of simuls for each SA value
nb_per_SA 
## [1] 100
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")
#To check 
# dim(sim)
# dim(sim_select_nongenetic)
# tapply(sim_select_nongenetic$index_SA_Gen, sim_select_nongenetic$SA_Gen, length)
# tapply(sim_select_nongenetic$index_SA_NonGen, sim_select_nongenetic$SA_NonGen, length)

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))


colors<-c("LMER" = "#DCA237", "Anova" = "#459B76")
#Plot
plot_nongenetic_binomial<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) + 
   geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   ggtitle("Binomial data") +
   theme_LO_sober +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_nongenetic_binomial

#Save plot
# name_plot<-paste0("Simul_powertest_nongenetic_binomial_",
#       nb_simul,"nbsimul_",nb_per_SA,
#       "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   plot_nongenetic_binomial, base_height = 10/cm(1), base_width = 14/cm(1), dpi = 1200)


#################################
### Common plot
#################################

power_plot_binomial<-cowplot::plot_grid(plot_nongenetic_binomial,plot_genetic_binomial,
                          labels=c("A", "B"), 
                          label_size = 15,
                          ncol =1, nrow = 2,
                     #    hjust = 0, vjust = 1,
                          scale = c(1,  1))
power_plot_binomial

#Save plot
# name_plot<-paste0("Simul_powertest_binomialdata_",
#       nb_simul,"nbsimul_",nb_per_SA,
#       "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   power_plot_binomial, base_height = 20/cm(1), base_width = 14/cm(1), dpi = 1200)
# 
# 

5.3 All plots

power_plot_all<-cowplot::plot_grid(plot_nongenetic_normal,plot_nongenetic_poisson,plot_nongenetic_binomial,
                                   plot_genetic_normal, plot_genetic_poisson,plot_genetic_binomial,
                          #labels=c("A", "B"), 
                          #label_size = 15,
                          ncol =3, nrow = 2,
                     #    hjust = 0, vjust = 1,
                          scale = c(1,  1))
power_plot_all

#Save plot
name_plot<-paste0("Simul_powertest_ALLDATA_",
      nb_simul,"nbsimul_",nb_per_SA,
      "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   power_plot_all, base_height = 20/cm(1), base_width = 42/cm(1), dpi = 1200)
# 

5.4 False positive rate analyses

### No genetic and no non genetic effects
nb_simul = 200
sim <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",  rho=0, rho_ng=0, 
              npop_per_fruit = 5, nfruit = 3, nhab = 3, nrep = 10, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
  
## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.01
mean(sim$IndicGen_aov)
## [1] 0.01
## Local adaptation (non genetic)
mean(sim$IndicNonGen)
## [1] 0.015
mean(sim$IndicNonGen_aov)
## [1] 0.015
### Genetic effects and no non-genetic effects
sim <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",rho=-0.1, rho_ng=0,
              npop_per_fruit = 5, nfruit = 3, nhab = 3, nrep = 10, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
tapply(sim$IndicGen, sim$SA_Gen, mean)
##        0.2        0.3        0.4        0.5        0.6        0.7        0.8 
## 0.00000000 0.16666667 0.06666667 0.45000000 0.17948718 0.23333333 0.26923077 
##        0.9          1        1.1        1.2        1.3 
## 0.39130435 0.36363636 0.50000000 1.00000000 1.00000000
tapply(sim$IndicGen_aov, sim$SA_Gen, mean)
##        0.2        0.3        0.4        0.5        0.6        0.7        0.8 
## 0.00000000 0.16666667 0.06666667 0.45000000 0.17948718 0.23333333 0.26923077 
##        0.9          1        1.1        1.2        1.3 
## 0.39130435 0.36363636 0.50000000 1.00000000 1.00000000
## Local adaptation (non genetic)
mean(sim$IndicNonGen)
## [1] 0.015
mean(sim$IndicNonGen_aov)
## [1] 0.015
### 
### Non genetic effects and no genetic effects
sim <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",rho=0, rho_ng=-0.1,
              npop_per_fruit = 5, nfruit = 3, nhab = 3, nrep = 10, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)

## Local adaptation (genetic)
mean(sim$IndicGen)
## [1] 0.07
mean(sim$IndicGen_aov)
## [1] 0.07
## Local adaptation (non genetic)
tapply(sim$IndicNonGen, sim$SA_NonGen, mean)
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.0000000 0.2500000 0.1000000 0.1875000 0.3061224 0.3783784 0.5000000 0.6000000 
##         1       1.1       1.2 
## 0.2500000 0.5000000 1.0000000
tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)
##       0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9 
## 0.0000000 0.2500000 0.2000000 0.1875000 0.3061224 0.3783784 0.5000000 0.6000000 
##         1       1.1       1.2 
## 0.2500000 0.5000000 1.0000000

6 UNBALANCED DESIGN

#Dataset

#unbalanced dataset
data_PERF <- read.table(file="/Users/laure/Documents/RESEARCH/z_PhD/POPNAT/OLD/DATA_COMPLET/Clean/DATACOMPLET_PERF_POPNAT2018.csv",sep=";",header=T)
data_PERF$Original_environment<-as.factor(data_PERF$Original_environment)
data_PERF <- data_PERF[data_PERF$Original_environment!="WT3",]
data_PERF <- data_PERF[data_PERF$Test_environment!="GF",]
data_PERF <- data_PERF[data_PERF$Test_environment!="Grape",] #mvrnom doesn't work with unbalanced matrix
data_PERF <- droplevels(data_PERF)

7 Normal distribution

7.1 Simulation

#Simul data 
nb_simul <- 10000
sim <- sapply(1:nb_simul, simul_fitnessdata_unbalanced, distrib = "normal", 
              unbalanced_dataset = data_PERF, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, rho = -0.05, rho_ng = -0.05, 
              sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)

7.2 Plot

####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5 
##   10  107  496 1127 1633 1919 1672 1307  868  464  234  106   37   19    1
#Save the same number of simuls for each SA value
nb_per_SA=100
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))


#Plot
plot_genetic_normal_unbalanced<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Normal data (unbalanced experimental design)") +
   theme_LO_sober +
   labs(color = "Methods", linetype = "Design") +
   scale_color_manual(values = colors) +
   scale_linetype_manual(values = lines)
plot_genetic_normal_unbalanced

# Balanced + unbabacaled
plot_genetic_normal_ALL <- plot_genetic_normal +
   geom_point(data = data_prop, 
              aes(y=prop_val_sign, x=val_SA_genet, color="LMER"), size = 2.5) + 
   geom_line(data = data_prop, 
             aes(y=prop_val_sign, x=val_SA_genet, 
                 color="LMER", linetype = "Unbalanced"), size = 1) + 
   geom_point(data = data_prop,
              aes(y=prop_val_sign_aov, x=val_SA_genet, 
                  color="Anova"), size = 2.5) + 
   geom_line(data = data_prop, 
             aes(y=prop_val_sign_aov, color="Anova", x=val_SA_genet, 
                 linetype = "Unbalanced"), size = 1) +
   labs(color = "Methods", linetype = "Design") +
   scale_linetype_manual(values = lines)
plot_genetic_normal_ALL

################################# 
### Plot non-genetic
#################################
#Number of simuls for each SA non genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.7 
##    4  111  493 1054 1759 1912 1695 1232  853  451  267   98   52   13    5    1
#Save the same number of simuls for each SA value
nb_per_SA 
## [1] 100
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))


#Plot
plot_nongenetic_normal_unbalanced<-ggplot2::ggplot(data = data_prop, 
                                                   aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Normal data (unbalanced experimental design)") +
   labs(color = "Methods") +
   scale_color_manual(values = colors)
plot_nongenetic_normal_unbalanced

plot_nongenetic_normal_ALL <- plot_genetic_normal +
   geom_point(data = data_prop, 
              aes(y=prop_val_sign, x=val_SA_nongenet, color="LMER"), size = 2.5) + 
   geom_line(data = data_prop, 
             aes(y=prop_val_sign, x=val_SA_nongenet, 
                 color="LMER", linetype = "Unbalanced"), size = 1) + 
   geom_point(data = data_prop,
              aes(y=prop_val_sign_aov, x=val_SA_nongenet, 
                  color="Anova"), size = 2.5) + 
   geom_line(data = data_prop, 
             aes(y=prop_val_sign_aov, color="Anova", x=val_SA_nongenet, 
                 linetype = "Unbalanced"), size = 1) +
   labs(color = "Methods", linetype = "Design") +
   scale_linetype_manual(values = lines)
plot_nongenetic_normal_ALL

8 Poisson distribution

8.1 Simulation

#Simul data 
#nb_simul <- 1000
sim <- sapply(1:nb_simul, simul_fitnessdata_unbalanced, distrib = "poisson", 
              unbalanced_dataset = data_PERF, sdpop = 0.5, 
              sdfruithab = 1, sdfruithab_ng = 1, rho = -0.05, rho_ng = -0.05, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)

8.2 Plot

####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9   2 2.1 
##   3  10  50 121 253 404 558 709 819 939 955 915 814 777 652 531 465 313 221 171 
## 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.1 
## 119  63  60  36  15  12  11   3   1
#Save the same number of simuls for each SA value
nb_per_SA=100
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))


#Plot
plot_genetic_poisson_unbalanced<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Poisson data (unbalanced experimental design)") +
   theme_LO_sober  + 
   labs(color = "Methods") +
   scale_color_manual(values = colors)
plot_genetic_poisson_unbalanced

plot_genetic_poisson_ALL <- plot_genetic_poisson +
   geom_point(data = data_prop, 
              aes(y=prop_val_sign, x=val_SA_genet, color="LMER"), size = 2.5) + 
   geom_line(data = data_prop, 
             aes(y=prop_val_sign, x=val_SA_genet, 
                 color="LMER", linetype = "Unbalanced"), size = 1) + 
   geom_point(data = data_prop,
              aes(y=prop_val_sign_aov, x=val_SA_genet, 
                  color="Anova"), size = 2.5) + 
   geom_line(data = data_prop, 
             aes(y=prop_val_sign_aov, color="Anova", x=val_SA_genet, 
                 linetype = "Unbalanced"), size = 1) +
   labs(color = "Methods", linetype = "Design") +
   scale_linetype_manual(values = lines)
plot_genetic_poisson_ALL

################################# 
### Plot non-genetic
#################################
#Number of simuls for each SA non genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)
## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9   2 2.1 
##   2  11  49 129 240 379 522 727 885 969 951 905 845 743 663 496 405 325 224 183 
## 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9   3 3.3 
## 127  78  53  35  26  15   5   3   4   1
#Save the same number of simuls for each SA value
nb_per_SA 
## [1] 100
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))


colors<-c("LMER" = "#DCA237", "Anova" = "#459B76")
#Plot
plot_nongenetic_poisson_unbalanced<-ggplot2::ggplot(data = data_prop, 
                                                   aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Poisson data (unbalanced experimental design)") +
   labs(color = "Methods") +
   scale_color_manual(values = colors)
plot_nongenetic_poisson_unbalanced

plot_nongenetic_poisson_ALL <- plot_nongenetic_poisson +
   geom_point(data = data_prop, 
              aes(y=prop_val_sign, x=val_SA_nongenet, color="LMER"), size = 2.5) + 
   geom_line(data = data_prop, 
             aes(y=prop_val_sign, x=val_SA_nongenet, 
                 color="LMER", linetype = "Unbalanced"), size = 1) + 
   geom_point(data = data_prop,
              aes(y=prop_val_sign_aov, x=val_SA_nongenet, 
                  color="Anova"), size = 2.5) + 
   geom_line(data = data_prop, 
             aes(y=prop_val_sign_aov, color="Anova", x=val_SA_nongenet, 
                 linetype = "Unbalanced"), size = 1) +
   labs(color = "Methods", linetype = "Design") +
   scale_linetype_manual(values = lines)
plot_nongenetic_poisson_ALL

9 Binomial distribution

9.1 Simulation

#Simul data 
#nb_simul <- 1000
sim <- sapply(1:nb_simul, simul_fitnessdata_unbalanced, distrib = "binomial", 
              unbalanced_dataset = data_PERF, sdpop = 0.5, 
              sdfruithab = 0.5, sdfruithab_ng = 0.5, 
              rho = -0.05, rho_ng = -0.05, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)

9.2 Plot

####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5 
##   10  107  496 1127 1633 1919 1672 1307  868  464  234  106   37   19    1
#Save the same number of simuls for each SA value
nb_per_SA=100
sim_select_genetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)), 
                 prop_val_sign = tapply(sim_select_genetic$IndicGen, 
                                        sim_select_genetic$SA_Gen, mean), 
                 prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov, 
                                        sim_select_genetic$SA_Gen, mean))


#Plot
plot_genetic_Binomial_unbalanced<-ggplot2::ggplot(data = data_prop, 
                                                  aes(x=val_SA_genet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant genetic local adaption")  +
   xlab("Intensity of genetic local adaptation (SA)")  +
   ggtitle("Binomial data (unbalanced experimental design)") +
   theme_LO_sober  + 
   labs(color = "Methods") +
   scale_color_manual(values = colors) 
plot_genetic_Binomial_unbalanced

plot_genetic_binomial_ALL <- plot_genetic_binomial +
   geom_point(data = data_prop, 
              aes(y=prop_val_sign, x=val_SA_genet, color="LMER"), size = 2.5) + 
   geom_line(data = data_prop, 
             aes(y=prop_val_sign, x=val_SA_genet, 
                 color="LMER", linetype = "Unbalanced"), size = 1) + 
   geom_point(data = data_prop,
              aes(y=prop_val_sign_aov, x=val_SA_genet, 
                  color="Anova"), size = 2.5) + 
   geom_line(data = data_prop, 
             aes(y=prop_val_sign_aov, color="Anova", x=val_SA_genet, 
                 linetype = "Unbalanced"), size = 1) +
   labs(color = "Methods", linetype = "Design") +
   scale_linetype_manual(values = lines)
plot_genetic_binomial_ALL

################################# 
### Plot non-genetic
#################################
#Number of simuls for each SA non genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)
##  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9    1  1.1  1.2  1.3  1.4  1.5  1.7 
##    4  111  493 1054 1759 1912 1695 1232  853  451  267   98   52   13    5    1
#Save the same number of simuls for each SA value
nb_per_SA 
## [1] 100
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim, 
                                          type_SA = "non-genetic",
                                          nb_simul_per_SA = nb_per_SA, 
                                          minimum_SA = "0")

#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)), 
                 prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen, 
                                        sim_select_nongenetic$SA_NonGen, mean), 
                 prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov, 
                                        sim_select_nongenetic$SA_NonGen, mean))


colors<-c("LMER" = "#DCA237", "Anova" = "#459B76")
#Plot
plot_nongenetic_Binomial_unbalanced<-ggplot2::ggplot(data = data_prop, 
                                                   aes(x=val_SA_nongenet)) + 
   geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") + 
   geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) + 
   geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") + 
   ylab("Proportion of simulations with\nsignificant non-genetic local adaption")  +
   xlab("Intensity of non-genetic local adaptation (SA)")  +
   theme_LO_sober + 
   ggtitle("Binomial data (unbalanced experimental design)") +
   labs(color = "Methods") +
   scale_color_manual(values = colors)
plot_nongenetic_Binomial_unbalanced

plot_nongenetic_binomial_ALL <- plot_nongenetic_binomial +
   geom_point(data = data_prop, 
              aes(y=prop_val_sign, x=val_SA_nongenet, color="LMER"), 
              size = 2.5) + 
   geom_line(data = data_prop, 
             aes(y=prop_val_sign, x=val_SA_nongenet, 
                 color="LMER", linetype = "Unbalanced"), size = 1) + 
   geom_point(data = data_prop,
              aes(y=prop_val_sign_aov, x=val_SA_nongenet, 
                  color="Anova"), size = 2.5) + 
   geom_line(data = data_prop, 
             aes(y=prop_val_sign_aov, color="Anova", x=val_SA_nongenet, 
                 linetype = "Unbalanced"), size = 1) +
   labs(color = "Methods", linetype = "Design") +
   scale_linetype_manual(values = lines)
plot_nongenetic_binomial_ALL

10 All plots unbalanced

power_plot_all_unbalanced<-cowplot::plot_grid(plot_nongenetic_normal_unbalanced,
                                   plot_nongenetic_poisson_unbalanced,
                                   plot_nongenetic_Binomial_unbalanced,
                                   plot_genetic_normal_unbalanced,
                                   plot_genetic_poisson_unbalanced,
                                   plot_genetic_Binomial_unbalanced,
                          #labels=c("A", "B"), 
                          #label_size = 15,
                          ncol =3, nrow = 2,
                     #    hjust = 0, vjust = 1,
                          scale = c(1,  1))
power_plot_all_unbalanced

#Save plot
name_plot<-paste0("Simul_powertest_ALLDATA_UNBALANCED",
      nb_simul,"nbsimul_",nb_per_SA,
      "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   power_plot_all_unbalanced, base_height = 20/cm(1),
#                   base_width = 42/cm(1), dpi = 1200)
# 

11 ALLPLOT

POWER_ALL<-cowplot::plot_grid(plot_nongenetic_normal,
                                          plot_nongenetic_poisson,
                                          plot_nongenetic_binomial,
                                          plot_genetic_normal,
                                   plot_genetic_poisson,
                                   plot_genetic_binomial,
                                   plot_nongenetic_normal_unbalanced,
                                   plot_nongenetic_poisson_unbalanced,
                                   plot_nongenetic_Binomial_unbalanced,
                                   plot_genetic_normal_unbalanced,
                                   plot_genetic_poisson_unbalanced,
                                   plot_genetic_Binomial_unbalanced,
                          #labels=c("A", "B"), 
                          #label_size = 15,
                          ncol =3, nrow = 4,
                     #    hjust = 0, vjust = 1,
                          scale = c(1,  1))
POWER_ALL

#Save plot
# name_plot<-paste0("Simul_powertest_ALLDATA_BALANCED_and_UNBALANCED",
#       nb_simul,"nbsimul_",nb_per_SA,
#       "nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
#                   POWER_ALL, base_height = 40/cm(1), 
#                   base_width = 42/cm(1), dpi = 1200)
# 



PLOT_ALL<-cowplot::plot_grid(plot_nongenetic_normal_ALL,
                                   plot_nongenetic_poisson_ALL,
                                   plot_nongenetic_binomial_ALL,
                                   plot_genetic_normal_ALL,
                                   plot_genetic_poisson_ALL,
                                   plot_genetic_binomial_ALL,
                          #labels=c("A", "B"), 
                          #label_size = 15,
                          ncol =3, nrow = 2,
                     #    hjust = 0, vjust = 1,
                          scale = c(1,  1))
PLOT_ALL

#Save plot
name_plot<-paste0("Simul_powertest_",
      nb_simul,"nbsimul_",nb_per_SA,
      "nbsimulperSA.pdf")
cowplot::save_plot(file = here::here("figures",name_plot ),
                  PLOT_ALL, base_height = 20/cm(1), 
                  base_width = 42/cm(1), dpi = 1200)